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ABSTRACT 

We solve for the time-dependent dynamics of axisymmetric, general relativistic magnetohy- 
drodynamic winds from rotating neutron stars. The mass loss rate as a function of latitude 
is obtained self-consistently as a solution to the MHD equations, subject to a finite thermal 
pressure at the stellar surface. We consider both monopole and dipole magnetic field geome- 
tries and we explore the parameter regime extending from low magnetization (low-o"), almost 
thermally-driven winds to high magnetization (high-cro), relativistic Poynting-flux dominated 
outflows (cr - /4npyc^ff x o-q/joo, P - vjc with ctq = w^O^/M, where u) is the ro- 
tation rate, O is the open magnetic flux, and M is the mass flux). We compute the angular 
momentum and rotational energy loss rates as a function of cry and compare with analytic 
expectations from the classical theory of pulsars and magnetized stellar winds. In the case of 
the monopole, our high-erg calculations asymptotically approach the analytic force-free limit. 
If we define the spindown rate in terms of the open magnetic flux, we similarly reproduce the 
spindown rate from recent force-free calculations of the aligned dipole. However, even for ctq 
as high as ~ 20, we find that the location of the F-type point {ry), which specifies the radius of 
the last closed field line in the equatorial plane, is not the radius of the light cylinder/?/, - cjo) 
(R = cylindrical radius), as has previously been assumed in most estimates and force-free cal- 
culations. Instead, although the Alfven radius at intermediate latitudes quickly approaches Rl 
as ctq exceeds unity, ry remains significantly less than Rl. In addition, ry is a weak function 
of (To, suggesting that high magnetizations may be required to quantitatively approach the 
force-free magnetospheric structure, with ry - Rl. Because ry < Rl, our calculated spindown 
rates thus exceed the classic "vacuum dipole" rate: equivalently, for a given spin-down rate, 
the corresponding dipole field is smaller than traditionally inferred. In addition, our results 
suggest a braking index generically less than 3. We discuss the implications of our results for 
models of rotation-powered pulsars and magnetars, both in their observed states and in their 
hypothesized rapidly rotating initial states. 

Key words: magnetohydrodynamics (MHD) — stars: winds, outflows — relativity — stars: 
neutron — methods: numerical 



1 INTRODUCTION 

Magnetically dominated winds from stars and accretion disks are 
central to the angular momentum evolution of these objects. Be- 
cause they can efficiently extract rotational energy — transforming 
stored gravitational binding energy into asymptotic wind kinetic 
energy — magnetic outflows ar e ubiquitous in pow ering a wide va- 
riety of astrophysical systems. Schatzman ( 1962) first introduced 
the key idea that a magnetic field anchored in a rotating star with 
a wind can enforce plasma co-rotation out to radii large compared 
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to the stellar radius, thus greatly increasing the angular momentum 
lost per unit mass. 

Schatzman's ideas were formalized in the steady non- 
relativistic flow model of IWeber & Davi j ^963), who assumed a 
pure monopole field geometry, and then by Mestel (1968ab) who 
assumed a more realistic dipole magnetic field configuration. For 
strong dipole fields, a closed zone forms at low latitudes and the 
mass loss is concentrated at high latitudes in regions where the field 
lines can be opened. If the extent and shape of the open field line 
region is known, then the physics of the magnetic wind is simi- 
lar to that in the monopole geometry: the flow emerges along the 
open flux tubes and its character is determined by passing through 
the slow magnetosonic (SM), Alfven (AL), and fast magnetosonic 
(FM) surfaces. If the thermal sound speed is small, the flow is 
driven primarily by magneto-centrifugal forces and the asymptotic 
outflow velocity is approximately va(Ra) ~ Ra<^, where R refers to 
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the cylindrical radius, oj is tlie stellar angular velocity, and Ra is the 
AL radius. Because Ra is the point at which B^(RA)/4n = pv^(Ra), 
the outflow kinetic energy is roughly equal to the magnetic energy. 
This remains true in the asymptotic wind. 

The theory of magnetic winds has been greatly extended and 
applied in its non-re lativistic form t o a wide variety of problems in 
stellar astrophysics. iGoslind h99d) and lAschwanden et alJ ilOOll) 
provide reasonably modem reviews in the solar context, where 
many of the recent developments in non-relativistic wind modeling 
have occurred. Multi-dimensional models have been motivated in 
part by the fact that toroidal magnetic fields can collimate the flow 
along the rotation axis via hoop stress into jets, a topic of much 
interest in t he modeling o f the bipolar outflows observed from pro- 
tostars (e.g.. Smith 1998). 

The discovery of rotation powered pulsars motivated the ex- 
tension of these ideas to relativistic flows and relativistically strong 
magnetic fields (B^IAn » pc^ at the source), in order to ac- 
count for the obse rved spin down of pulsars. Michel ( 1969) and 
iGoldreich & JulianI |1970) made the first extensions of the Weber 
& Davis model to relativistic outflows, recovering the same se- 
quence of critical points in the flow. As reviewed in ^ these stud- 
ies yielded estimates of the spindown torque of 



k~ 1. 



(1) 



The dipole moment of the star, yu, was assumed to be related to the 
monopole moment m appearing in the MHD theory via m = fi/R^, 
a relation equivalent to assuming that the AL radius Ra was equal 
to the radius of the light cylinder Ri = c/ai. Since these monopole 
models only treated the flow in the rotational equator, the numerical 
coefficient ~ 1 in equation Q could not be determined accurately. 

These simplified theoretical models revealed important dif- 
ferences between relativistic and non-relativistic winds. First, in- 
stead of reaching approximate energy equipartition between flow 
kinetic energy and magnetic energy, the asymptotic flow remains 
strongly magnetized. The asymptotic Lorentz factor is given by 
7co ~ 0"^^, where ctq = 9^ci//Mc = B^^jAnpc^, O is the mag- 
netic flux, and M is the mass loss rate. Thus, ctq. For a 
highly relativistic outflow (ctq » 1) the asymptotic magnetiza- 
tion (Too - ctq/Joc ~ o"q'^^ s> la second important difference 
is that in relativistic magnetized flows the electric force cannot be 
neglected' . Although it is absent by symmetry in the Michel and 
Goldreich & Julian models, the electric force almost exactly can- 
cels the focusing hoop stress in multi-dimensional monopole mod- 
els, thus undermining these flows' usefulness for the understanding 
of relativistic col limated jets. This is a probl em generic to all rela- 
tivistic outflows jLvubarskv & Eichlej200ll) not focused by some 
external medium (e.g., a disk or external channel). 



' In relativistic MHD the electric field arises because of motion of the con- 
ducting fluid across the magnetic field. Equivalently, the fluid moves with 
the single particle ExB drift velocity. Therfore, one can eliminate the elec- 
tric field by puting oneself in the local fluid rest frame and describing the 
electromagnetic stress as being solely due to the comoving magnetic field. 
Occasionally, debate occurs about which is the "correct" frame in which to 
describe the forces, since workers used to non-relativistic MHD are some- 
times uncomfortable with the explicit appearance of the electric field. When 
the fluid dynamics is expressed in covariant form (the electrodynamics is al- 
ready covariant), as we have done in this paper, the forces are well defined 
and unambiguously described in every reference frame. The choice is dic- 
tated only by convenience in describing the physics, or in actually carrying 
out the calculations. 



Much effort has gone into relaxing the simplifications of these 
early models, and in particular, on understanding what is required 
for an ideal MHD flow to have y^o — > o"o- Observations and mod- 
els of pulsar wind nebulae suggest that at the termination shock 
of young rotation-powered pulsars' outflows the magnetic energy 
has been almost completely converted into flow kinetic energy (the 
"tr problem") and the flow 4-velocity has reached y a; cto (the "y 
problem"), in contradiction with the predictions of the one dimen- 
sional monopole treatments. If the magnetic surfaces retain an al- 
most m onopolar shape the ac celeration of the flow is only loga- 
rithmic iBegelman & L|1994. However various asymp totic solu- 
tions teegelman & L ll994lHevv aerts & Norman"2003' and refer- 
ences therein) and similarity solutions ( Vlahakis & Konisl 200^, 
Vlahakis 2004) of parts of the full axisymmetric flow problem sug- 
gest that if the poloidal magnetic field is sufficiently non-monopolar 
beyond the FM point, the largely radial current that supports the 
toroidal magnetic field, might cross field lines toward the equa- 
tor, causing a transfer of electromagnetic energy to flow kinetic 
energy. However, both a numeri cal solution in axisymmetric non- 
relativistic MHD iSakuraill985h and a perturbative calculation of 
relativi stically magnetized MHD flow from the exact force-free so- 
lution jBeskin et alJll993) yield a poloidal magnetic field which 
hews closely to the exact monopolar form. Thus, the monopole 
puzzle of asymptotic winds with magnetic energy that is never con- 
verted to flow kinetic energy persists, suggesting that ideal MHD 
expansion of the wind contradicts the observations. 

The monopole geometry for the poloidal field has long been 
recognized as a singular case. A dipole magnetic field aligned with 
the stellar rotation axis presents the simplest realistic alternative 
field geometry. Because the open field, where outflow occurs, has a 
nozzle shape that expands faster than at distances comparable to 
the radius of the last closed field line (ry), there is a possibility of 
faster acceleration and larger conversion of magnetic energy into 
kinetic energy than occurs in the monopole geometry. However, 
complete solutions in axisymmetric steady flow have only been 
possible for the monopole; the change of topology between closed 
and open field lines required in the dipole case has so far escaped 
solution of the Grad-Shafranov equation that describes the mag- 
netic structure in MHD flows with inertia and pressure included. 

The location of the AL radius in the outflow, and its relation to 
the maximum equatorial radius of the last closed field line ry, where 
typically a Y-type neutral point occurs in the magnetic field, is a 
problem of equal significance. As outlined in ^ if < R^, and if 
ry /Rl is an appropriate function of <i>, one might be able to account 
for the observed fact that in rotation powered pulsars Ty oc a)", with 
n < 3. The strong inferred magnetization of pulsars has always sug- 
gested that the outflow structure should be treated as a problem in 
force-free relativistic MHD, with inertial and pressure forces com- 
pletely neglected ( Goldreich & Julian 1969). Solution of the force- 
free Grad-Shafranov equation for the case o f an aligned dipole 
jScharlemarin & WagoneJl97llMichelll973h resisted solution un- 
til recently <Contopoulos et alj|l999l and lGruzinovlbOOSl assumed 
strictly force-free conditions evervwhere. whereas ^Good win et al] 
2004 included pressure effects near the Y-point). Gruzinov found 
that the torque is indeed given by equation Q, with k = 1 ± 0.1, 
while Contopoulos et al. found k = 1.85. All of these authors as- 
sumed ry = Rl, as is physically plausible. 

T, 



However, subsequent work by Timokhin ( 2005) has found that 
force-free MHD allows a family of solutions with ry of the last 
field line within Ri, contrary to naive expectations. He found that 
k x QAliRtlry)^ , which, when ry = Rl, is close to Michel's force- 
free monopole result. Importantly, these solutions indicate that be- 
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yond a few Light Cylinder radii, tlie poloidal magnetic field as- 
sumes a monopolar form, which suggests that acceleration in the 
asymptotic wind never reaches equipartition energies (the cr prob- 
lem). However, since inertial forces are important beyond the FM 
surface (r « ct'J^Rl) the force- free approximation breaks down 
jBeskin et al]ll998t lAronsll2004 and no strong conclusions may 
be drawn. 

The ambiguities of the steady force-free models, and the dif- 
ficulty of solving the magnetospheric structure in MHD with iner- 
tia and pressure included, suggest that an evolutionary approach 
to the problem would be useful. A time-dependent solution can 
then be sought for specified boundary conditions at the stellar sur- 
face (such as pressure or injection velocity). The system can then 
relax to a self-consistent steady state (if one exists) and this ap- 
proach allows one to find the last closed flux surface unambigu- 
ously, for specified injection conditions. Furthermore, such a treat- 
ment allows for the possibility of intrinsically time-dependence 
(including either limit-cycle or chaotic behavior) of the magne- 
tosphere. The time dependence of individual radio pulses from 
rotation-powered pulsars, showing systematic drifting through the 
pulse window for many objects near the pulsar death line, and 
chaotic rotation phase over wider regions of the P - P diagram 
(Rankin 1986, Deshoandhe & Rankin 1999) have timescales con- 
sistent with global magnetospheric variability causin g changes in 
the polar current system underlying pulsar emission iArondl 19811: 
[Wri aht 2001). I ndeed, th e torque noise exhibited by many pulsars 
jCor des & Helfand 1980) may also owe its origin to instabilities of 
the magnetospheric current system I Arons 1981. ; , Cheng 1987a b,). 

In this paper we carry out time-dependent RMHD modeling of 
both highly magnetized Poynting-flux dominated winds (ctq » 1) 
in which Ra ~ Rl, as well as magnetized neutron star winds in 
which Ra is significantly less than Ri^ (ctq <k 1). 

Large ctq MHD models draw their motivation from studies 
of rotation-powered pulsars and of magnetars (Soft Gamma Ray 
Repeaters and Anomalous X-ray Pulsars). More broadly, analo- 
gous problems appear in Poynting-flux dominated models of jets 
from black holes and neutron stars. Low-cto (but still magnetized) 
neutron star winds are of interest primarily in understanding the 
physics of very young neutron stars. In the seconds after collapse 
and explosion, the neutron star is hot (surface temperatures of ~2-5 
MeV) and extended. This proto-neutron star cools and contracts 
on its Kelvin-Helmholtz timescale (tkh ~ 10 - 100 s), radiating 
its gravitational binding energy (~10^^ ergs) in neutrinos of all 
species furrows & Lattimer 1986, Pons et al. 1999). The cooling 
epoch is accompanied by a thermal wind, driven by neutrino en- 
ergy deposition (primarily Vi,n — > pe" and v^p — » ne^), which 
emerges into the post-supemova-shock ejecta (e.g . iDuncan et alJ 
ll986llOian & WoosIevll9MlThompson et all200ll) . 

A second or two after birth, the thermal pressure at the edge 
of the protoneutron star surface, where the exponential atmosphere 
joins the wind (r^), is of order ~ 10^^ ergs cm"^ and decreases 
sharply as the neutrino luminosity (L,,) decreases on a timescale 
Tkh, as the star cools and deleptonizes. The thermal pressure at 
the stellar surface is set by Ly. If the protoneutron star has a sur- 
face magnetic field of strength So, then at some point during the 
cooling epoch the magnetic energy density will exceed the ther- 
mal pressure. For fixed So, the wind region becomes increasingly 
magnetically dominated as Ly decreases. For larger So the mag- 
netic field dominates at earlier times. For magnetar- strength surface 
fields (So ~ lO'** - 10'^ G) the magnetic field dominat es the wind 
dynam ics from just a few seconds after the supernova iThompsorl 



Magnetars are thought to be bom with millisecond rotation pe- 
riods jDuncan & Thompsonlll993 . lThomDson & DuncaDll993ll . in 

which case the combination of rapid rotation and strong magnetic 
fields makes the proto-neutron star wind magneto-centrifugally 
driven. Because the rotational energy of a millisecond magnetar 
is very large (~2x 10^^ ergs) relative to the energy of the supernova 
explosion (~10" ergs) and because the spindown timescale Tj ~ 
ci)/cij ~ {2I5){MI M){rNs I rA?' (where r^s is the neutron star radius) 
can be short for large va, these proto-magnetar winds have been 
consid ered as a mechanism f or producing hyper-energetic super- 
novae iThompson et all2004) . Because the wind becomes increas- 
ingly magnetically-dominated and the flow eventually becomes 
Poynting-flux dominated as the neutrino luminosity abates, the 
outflow must become relativistic. For this reason, proto-magnetar 
winds have also been considered as a pos sible centra l engine for 
long-duration gamma ray bursts (GRBs) jysovlll992L iThompsor] 
[1994', 'Wheeler et al.''2000', (Thompson et alj|2004 '). They may also 
be the source of ultra-high energy cosmic rays iBlasi et al. 200^ 
IXrons 2003). The time dependent RMHD models of neutron star 
winds calculated in this paper provide a significantly improved 
understanding of proto-magnetar winds, and their possible role in 
hyper-energetic supernovae and GRBs. 



1.1 This Paper 

With these goals in mind, in ^we first outline some order-of- 
magnitude estimates for the spindown of neutron stars in both the 
low- and high-cro limits. In ^& ^ we describe the details of our 
numerical model. In ^we present our numerical results for the ID 
monopole ( ^5.1> . the 2D monopole (' ^5.2> . and the aligned dipole 
( 35. 3t . In each case, we present a set of models for both low- and 
high-iTo winds. Finally, in ^we present a discussion of our results 
and speculate on their implications for the spindown of rotation- 
powered pulsars and very young, rapidly rotating magnetars. 



2 PHYSICAL MODEL 

As a guide to the numerical models, we describe here several sim- 
ple order of magnitude estimates of the properties of both high- and 
low-tTo winds from neutron stars. 



2.1 Poynting Flux Dominated Spindown 

Consider a star with a magnetic dipole moment ji. For simplic- 
ity, assume /j || uj. Suppose there is an outflow of plasma along 
open field lines which connect to the star in a polar cap, with 
the magnetic flux of the open field lines being O^. The expected 
poloidal magnetic structure is shown in Figure Q In the closed 
zone, plasma co-rotates, and the toroidal currents, composed of 
co-rotating charge density and pressure and inertial drifts across 
the magnetic field, cause the distortions from the "vacuum dipole" 
field, which are of importance at radii comparable to ry. Assume 
the AL radius Ra is comparable to ry. The Poynting flux is 5 = 
(c/4n)E X S, whose radial component is, with the poloidal electric 
field E = -(a) x r) x Bp, 5, = -(cursme/c)Bp. With Bp(rA) ~ 
Bd.poieirA) and S/r^ sin0) « -S„(i?^), S ,(Ra) ~ (wRA/47T)fi^ /R% 
where subscripts p and <p denote the poloidal and toroidal compo- 
nents, respectively. Then the EM spindown torque is approximately 
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with decreasing spin down power. One way to state this is to 
simply assert that Ra = ry < Rl and that ry/R^ oc ca" (e.g., 
''y/Rl = ftjfA/c /c')°'. lAronJl983h . This implies a change in the po- 
lar cap size from r^sirNS IRlY^^ to the larger value (larger Oo) of 
rNsirNS /Rl)^' "^'~- Using this expression, the braking index data 
require 1/6 > a > 1/30, with the largest value for the Crab pulsar 
and the smallest for the 407 ms pulsar Jl 119-6157. 

Assuming a > is equivalent to the last closed field line of 
the dipole having equatorial radius ry less than R^ ~ Rl- Accord- 
ing to our estimate J2}, the electromagnetic torque depends on the 
field strength at Ra, which is noticeably larger than that estimated 
by using a pure dipole filed, since the poloidal field becomes pro- 
gressively more monopolar for r > ry. Thus one obtains a better 
estimate of the torque by using a simplified model of the poloidal 
magnetic field which has the correct asymptotic form shown by the 
force-free aligned rotator models - dipolar at r <k ry and monopolar 
at r » ry. Our RMHD results have the same asymptotic behavior. 
Thus, with Bp = (ulr^) + K{p.lryr^), the same order of magnitude 
argument that led to J2j yields 



3 

, >^ hi' 



(3) 



Figure 1. Force-free poloidal magnetic field lines from a magnetized star 
with dipole axis aligned with the rotation axis. The distances are sca led in 
units of the radius of the Light Cylinder, Rlc- From .Gruzino\l l2005t) . 



Anrl 



(2/3)S,iRA) 



3~ 



(2) 



the evaluation of the constant k in equation Q to be equal to 2/3 
is the exact re sult for the EM torque on the force-free monopole 
jMiche]|l973ll . which has no closed zone or Y point. Since the open 
magnetic flux is <bg x (ij,ci)/c)(RL/ry), the torque can also be written 
as Til « {2/3mlca/c)(ry/RLf. 

If mass loading (inertial forces) and pressure are negligi- 
ble, the long standing expecta tion has been ry = Ra = Rl = 
c/a> iGoldreich & Juliaall969f) . and MHD spindown torques then 
should have braking index n = (bUj/ar = 3, if p and the stellar 
moment of inertia are constant. Two numerical solutions for steady 
flow from the force-free a ligned dipole rotator have found k = 1.85 
Conto goulos^ alJll99St) and, at higher resolution, k = 1 ± 0. 1 
Gruzinovll2005h . rather than 2/3. The equality ry = Rl = Ra was 
assumed in both of these calculations. 

A long standing empirical puzzle has been that in the four 
observations of braking indices not requiring major corrections 
for glitc hes in the timing , the braking index lies between 2.5 
and 2.9 jLvne et alJll993L iKaspi et alJll994 lOeeter et alJll999l 
ICamilo e^ L200?! TJvhtgston^ ^r*2^ reduction of the 

braking index for fixed yj and ai, in comparison to our simple es- 
timate, is tantamount to ry < Rl. That is, the closed zone ends 
within the Light Cylinder^, and, as the star ages (ai decreases), 
ry/RA also decreases — the magnetosphere becomes more open 



An alternate possibility is an increasing magnetic moment 
iBlandford & Romani Il988li . Still another option is evolution of the 
angle between cj and //, in the still unassessed dependence of the elec- 
tromagnetic torque on the oblique rotator, in either force-free or RMHD 
models. 



r _ "/I 

/ = k:—- 

ry 



(4) 



Assuming the magnetic moment and the stellar moment of inertia 
(and ; = /.(co,fi) are constant, and that Ra = Rl, a correspondence 
assumed in most force free models and also found in the RMHD 
results we report below, one readily finds 



Q.CL) 



3 + 2 



31n(l +/) 
dlno) 



(5) 



Thus a braking index less than 3 requires RaI^y to depend on oj 
(more generally, to depend on time.) If the time dependence of 
rylRA enters solely through dependence on uj, n < 3 requires 
RA/ry = RL/ry to increase as the star spins down - the magneto- 
sphere becomes more open with time, and the magnetic field at the 
light cylinder to remain larger, than is expected in the traditional 
model. Such behavior requires transformation of closed field lines 
to open, which can occur if magnetic dissipation at and near the 
Y-point allows reconnection to enable this transformation. 

Our RMHD numerical results presented in ^5. 3 1 show that al- 
though Ra — » Rl as soon as ctq exceeds unity, ry /Rl remains sub- 
stantially less than unity for <Tq as large as ~20. We also find that, 
for (To of order a few, the ratio ry//?/. decreases as a> decreases - 
the braking index in our models is less than 3. This suggests that 
seemingly small inertial and pressure forces can have a large effect 
on the magnetospheric structure and, in turn, the magnitude of the 
spindown torque and the br aking index. 

The work of Mestel & Spruill fl987i) may provide an explana- 
tion for our numerical results. In their isothermal, non-relativistic 
analysis of the magnetohydrostatic equilibrium of the closed zone, 
they find that Ry/rj^g depends crucially on the ratios Vll2c\. and 
a/ry/V^, where Vg is the escape velocity from the stellar surface 
(see their §2 & eq. 8). For ry a few times rj^g one finds an approxi- 
mate implicit equation for ry: 

ry\' ^ (Bl/U) 
rNS 



Pns c\ 



■ exp 



1/3 



2c\ \ VI 



(6) 



One sees that for fixed isothermal sound speed cj, if (wry/V^) is 
greater than unity, then it becomes exponentially harder to increase 
RYlrNS by increasing B\Ipns. 

However, since electromagnetic stresses alone can lead to Y 
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point formation, as is clear from the solutions of the force-free 
Grad-Shafranov equation, one should treat these estimates of Ry 
as a function of plasma stress alone with caution, when applied to 
the magnetically-dominated regime. Simplified models along the 
lines pioneered by Goodwin et al. (2004) may be helpful, but at 
present, the best results are the simulations th emselves (^5.3> . We 
defer a relativistic generalization of .Mestel & Spruij il987i) with 
a polytropic (p oc p^) equation of state (more appropriate to our 
simulations) to a future paper. 

2.2 Spindown by Magnetized Mass-Loaded Winds 

Stresses due to mass loading become significant in the thermally 
driven winds of newly born neutron stars. These stars may have 
strongly mass-loaded winds which have large in comparison to 
rjv5 , but significantly smaller than Rl ( Thompson et al. 2004). Con- 
ditions for such winds in the presence of thermally-driven mass 
loss generally obtain when the isothermal sound speed cj- is smaller 
than, or of order, the Alfven speed, va 

The torque on the star by the magnetized wind is easily es- 
timated for a strictly monopolar field geometry (Bp oc r^-). At 
we expect the poloidal magnetic energy density to be of order the 
radial kinetic energy density, B^Jin ^ pvJ/2, and if Cj is much 
smaller than v^, then we further expect that v^irA) should be of or- 
der v^(rA). A simple estimate for v^(rA), which again requires that 
Va '» ct, is ~ vaCi)- We combine the above ingredients to obtain 
an expression for the Alfven point: 

rA = Bl'V>^(Mcr'\ (7) 

where we have assumed that M = Anr^p\\ and So is the magnetic 
field strength at the stellar surface. The time evolution of the star's 
spin period is then governed by the equation 

j = -Mrloj = -B'J\'-II{MojYI\ (8) 

The rotational energy loss rate is £ = Iojlo oc M''^. In this limit 
(va » cj) the sonic point is approximately 

where we have scaled the spin period P to 1 millisecond and rsdi 
is the Schwarzschild radius. In the Cj- <k va limit, the radial veloc- 
ity r eaches its asymptotic value o f Voo ~ (3/2)va at the FM point 
fe.g. lBelcher & MacGregoJl97d) . 

In the regime we are interested in this paper we always have 
a supersonic outflows at the AL surface. In this case, according 
for the standar papa metrization of MHD winds (Sakurai 19851 
IPaigne & DrenMiSir>.200 2). all solution should be considered cen- 
trifugally or marginally centrifugally driven. 

A similar set of estimates for a dipole magnetic field is consid- 
erably more complicated, particularly since, as in ^2.11 the position 
of the Y point is not known a priori. In addition, since the areal 
fiinction along open flux tubes adjacent to the closed zone devi- 
ates significantly from radial, a full numerical solution is required 
to address spindown in this context 05.3> . However, as a rough 
guide in understanding the expected differences between monopole 
and dipole spindown it is sufficient to imagine the scalings for the 

^ For very high mass loss rates and/or high thermal pressures, pc~. s> 
B^/Sn. In this limit, the magnetic field is not important in accelerating mat- 
ter off of the stellar surface, the wind is driven thermally, and spindown is 
negligible unless the star rotates at a significant fraction of breakup. 



dipole as essentially those for the monopole, but with the surface 
magnetic field strength normalized to just the open magnetic flux. 



2.3 Parameterization 

These estimates and those of 82.11 reveal the parameters which 
specify the physical regimes of relevance to our models of rotating 
magnetospheres. Mass loss is thermally and centrifugally driven in 
these models, depending upon the ratio of pressure and centrifugal 
forces to the gravitational force, parameterized at the injection sur- 
face (r,„) by (cl/V^),-=ri„ and by (fir,,,/ V^,);^^, . All of the models 
considered in this paper have the first ofthese parameters between 
0.01 and 0.1, while the second is between 0.05 and 0.3. The val- 
ues adopted can be derived from the paramenter shown in Tables 
1 ( 84. U . For all of the 2D monopole and dipole models, magnetic 
pressure dominates gas pressure, as expressed by B^/4ncl > 1. 
This is true for most of the ID monopole models also. Again, these 
parameters are listed in Table 1. In all cases, the thermal energy 
density is smaller than the rest mass density (p < pc'). Thus, pres- 
sure forces do not lead to relativistic motion. The values of the ra- 
tios between the characteristic speeds at the base of the wind, for 
all our simulations, are provided in Appendix A. 

The distinction between pressure deriven and centrifugally 
driven wi nd, can be also done based o n the cond itions at the A L 
surface. In lDaigne & DrenkhahrJ<2002h (foUowing lSakurail < 1 985l) ) 
the distincion is based on the value of the ratio Tplp{Q.rY at the AL 
surfaces (F is the adiabatic coefficient). In all our cases such ratio 
is less that 0.1. 

The most significant parameter is Michel's magnetization pa- 
rameter, o"o, defined just after expression Q. When the magnetic 
energy density exceeds the rest mass density (cto > 1), magnetic 
pressure can accelerate the flow to relativistic velocities. This pa- 
rameter is listed for all the models, as our major goal is to span the 
regimes from highly mass loaded, pressure driven, nonrelativistic 
outflow (cTo <K 1) to lightly mass loaded, magnetically driven rel- 
ativistic outflow (cTo » 1) in both monopole and dipole geometry. 
The values of o-q for the various models are given in the tables in 
§ a5.1ll5.2l andl531 

We do not consider outflows driven by relativistically high 
temperature {p > pc^), a regime more relevant to fireball models 
of Gamma Ray Bursts. 



3 MATHEMATICAL FORMULATION 

The laws of mass and momentum-energy conservation, to- 

S ;ether with Maxwell e quations i n ge n eral relativity, are 
Landau & Lifshit3 ll97lL IWeinberd Il97l iMisner et alJ Il973l 
lAnile^l989^ : 

Vy(pU'') = 0, (10) 

v,(r'"') = o, (11) 

dpFy, + d,F^^ + d,iFpy = 0, (12) 
V^(F'") = -r, (13) 

where p is the proper rest mass density, m" and y are the four- 
velocity and the four-current density, and F'"' is the Faraday tensor 
of the electromagnetic field. The momentum-energy tensor is given 
by: 

r"" = phifu" + g^"p + F'"'F^^ - g^'F^'F^/A, (14) 
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where g*"" is the metric, p is the gas pressure, and we have chosen 
a system of units in which c = 1 . In the case of F-law equation of 
state for a perfect gas the specific enthalpy is = \+T l(r-\)pl p.\n 
order to close the system, the current density y must be specified in 
terms of the other known quantities, through an additional equation, 
Ohm's law. In the MHD approximation Ohm's law becomes the 
condition that the net electric field in the fluid frame must vanish, 
which in covariant form reads: 



F„„k'' = 0. 



(15) 



With this approximation, Eqs. <1U> - <13> can be rewritten in term of 
proper density, pressure, four-velocity and magnetic field, reduc- 
ing to a system of 8 equations for 8 variables, plus the solenoidal 
condition on the magnetic field, V • B = 0. 

Although we consider winds from neutron stars with centrifu- 
gal forces large enough to aff^ect the mass loss, we only consider 
rotation rates slow enough to allow us to neglect rotational modifi- 
cations of the metric. Therefore, employ the Schwarzschild instead 
of the Kerr metric. The use of a diagonal metric allows one to sim- 
plify the equations and to implement them easily in any code for 
special relativistic MHD, as shown bv lKoide et alJiT99^ . 

In Boyer-Lindquist coordinates {t,rfl,^) the diagonal elements 
of the metric are: 



-g 



=(l-2GM/r): 



(16) 



where M is the mass of the central object (1.44 Mg). 

In axisymmetry and steady state, equations <10> - J13> admit 
6 integrals of motion along stream lines (flux tubes ) (Camenzind 
1 1 986 iflCamenzin J 1 986bllCamenzin J 1 98llDaigne' & Drenkhahn 



T = apyVpA, 


(17) 


O = B„A, 


(18) 


n = a{v^ - B^/BpVp)/R, 


(19) 


£ = R{hyv, - a<S>B^/ri 


(20) 


•H = a{hy - Q.<S>RB^/r), 


(21) 


S = p/p\ 


(22) 



where subscript p and (p indicate poloidal and azimuthal compo- 
nents, respectively, y is the Lorentz factor, A is the area of the flux 
tube (A = in the case of radial outflow), and R = rsind the 
cylindrical radius. Again note that we have chosen units in which 
c = 1. 

If the value of the six integrals is known on a stream line, 
then equations <17> - <22t can be solved for the value of primitive 
quantities like density, velocity, and pressure. In general, quanti- 
ties like n,S, and O are assumed to be known and given by the 
physical conditions at the surface of the central object, while the 
values of the remaining integrals are derived by requiring the solu- 
tion to pass smoothly the slow-magnetosonic (SM), Alfvenic (AL) 
and fast-magnetosonic (FM) points. In 2D, one also requires an ad- 
ditional equation for the area of the flux tube, and this is provided 
by requiring equilibrium across streamlines. 



particularly simple and efficient, since solvers based on character- 
istic waves are avoided in favor of central-type component-wise 
techniques (HLL solver based only on the FM speed). In the axi- 
symmetric 2D approximation the equation for the evolution of B^ 
can be written in conservative form and only one component of 
the vector potential, A^, is used in the evolution of the poloidal 
magnetic field. Moreover we replace the energy conservation equa- 
tion with the constant entropy condition, S = const. Of course 
this condition cannot be satisfied during the evolution when shocks 
form in the flow. However, if the wind evolves toward steady state, 
shocks propagate outside the computational domain, and the con- 
dition S = const can be satisfied. There are various numerical rea- 
sons for our choice of 5 = const. Most importantly, enforcing con- 
stant entropy significantly enhances code stability. For example, it 
is known that in strongly magnetized flows or in the supersonic 
regime, in deriving the thermal pressure from the conserved quan- 
tities, small errors can lead to unphysical states. In non-relativistic 
MHD these unphysical states correspond to solutions with negative 
pressure. A common fix is to set a minimum pressure "floor" that 
allows the computation to proceed. However, in relativistic MHD 
(RMHD) it is possible that no state can be found (not even one with 
negative pressure) and the computation stops. This usually happens 
for high Lorentz factor y ~ 10- 100, depending on the grid-flow ge- 
ometry, or in the case of high magnetization, B^/(ph) > 100. When 
the magnetization at the Light Cylinder is high, close to the central 
object it will be above this stability threshold, which prevents us 
from using the full equation of energy. The use of a constant en- 
tropy condition also greatly simplifies the derivation of primitive 
variables from the conserved quantities, thus increasing the effi- 
ciency of the code. 

Lastly, by specifying the value of the entropy we remove en- 
tropy waves from the system. Entropy waves travel at the advection 
speed (and are dissipated in an advection time). If the sound speed 
is much smaller than c, the advection speed Vr, close to the surface 
of the central object, is much smaller than the speed of light and the 
advection time can be extremely long. By removing entropy waves, 
perturbations are dissipated at the SM speed, which in our regime 
is of order of 0.1 c. 

Simulations were performed on a logarithmic spherical grid 
with 200 points per decade in the radial direction, and a uniform 
resolution on the 9 direction with 100 grid points between the pole 
and the equator (CFL factor equal to 0.4). Higher resolutions were 
used in a few cases, to test convergence and accuracy. In order to 
model the heating and cooling processes using an ideal gas equa- 
tion of state we adopt an adiabatic coefficient F = 1.1 (almost 
isothermal wind) , which is reasonably re presentative of the wind 
solution by, e.g.. ioian & Wooslevj il996h . It can be shown that in 
order to have a transonic outflow the thermal pressure cannot be too 
high (above a critical value depending on F the sonic point moves 
inside the surface of the star) n or too lo w (the Bernoulli integral H 
must be positive), as shown bv lKoide'et al. (J 999) in the hydrody- 
namical case. The available parameter space increases for smaller 
F. The value we choose allows us to investigate easily cases with 
p/p ~ 0.05, especially in the dipole case where the strong flux di- 
vergence at the base is more important. 



4 NUMERICAL METHOD 

Equations <10> - <13> are solved using the shock-capturi ng code for 
relativistic MHD developed bv lDel Zarma et all i2003h . The code 
has been modified to solve the equations in the Sc hwarzschild met- 
ric following the recipes by Koide et al] ^99^. The scheme is 



4.1 Initial and boundary conditions, and Model Parameters 

Simulations in both ID and 2D were initialized by using a hydro- 
dynamical ID radial solution obtained on a much finer grid (the 
relativistic extension of the Parker solution jParkeJl95 8l). and pro- 
jecting it on the initial magnetic field lines. In the monopole cases. 
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initial poloidal magnetic field lines are assumed to be strictly ra- 
dial while for the dipole cases we adopt the solution for the vac- 
uum dipole in the Schwarzschild m etric (i.e.^uslimov & TsygaJ 
ll992llWasserman & Shapiro" 1983\ Density and pressure were in- 
terpolated from the hydrodynamical solution, v,- was derived by 
projecting the radial velocity on the magnetic field lines, and we 
set vg = 0. We also impose co-rotation in the inner region = 
mm(ClR/a, 0.6c), in order to avoid sharp temporal transients in the 
vicinity of the inner boundary. 

Standard reflection conditions are imposed on the axis, and 
symmetric conditions are imposed on the equator. At the outer ra- 
dial boundary we apply standard 0-th order extrapolation for all 
the variables. Initial condition are chosen in order to guarantee that 
during the evolution the FM surface is inside the computational 
domain so that no information is propagated back from the outer 
boundary. 

Unfortunately, as we will discuss in the following section, in 
the 2D case such a constraint can not be satisfied close to the axis, 
unless one uses an excessively large computational domain. In a 
few cases, using larger grids that allow the FM surface to be inside 
the computational domain, we find that the results do not change 
appreciably except along the axis itself. That is, the global solution 
at all but the highest latitudes is not significantly affected by the 
fact that the FM surface is outside the computational domain very 
near the polar axis. 

Particular care has to be taken for the inner boundary condi- 
tions. As pointed out in ^ we are here interested in the transi- 
tion from mass loaded (tro <c 1) to high-cro winds, and our in- 
jection conditions are tuned to be as close as possible to the neu- 
trino driven proto-neutron star case. We chose the inner radius r,„ 
to be located at II km (1.1 radii of the neutron star r^s ), which 
corresponds to the outer ed ge of the exponential a tmosphere for a 
thermally-driven wind (e.g. iThompson et all200lh . The modeling 
of such a steep atmosphere requires very high resolution in order to 
avoid numerical diffusion, and the problem becomes prohibitive in 
terms of computational time in 2D. At the inner boundary the flow 
speed is smaller than the SM speed, implying that all wave modes 
can have incoming and outgoing characteristics. This constrains the 
number of physical quantities that can be specified. Density and 
pressure at the inner radius are set to be p/p ~ 0.04, thus fixing 
the entropy for the overall wind. The radial velocity is derived us- 
ing linear extrapolation. We also fixed the value of n/Q-(r,„), typ- 
ically at nr,„/Qr(r,„) = 0.2, corresponding to a millisecond period 
(in the neutron star proper frame, see also Eq. <24> '). which implies 

The frozen-in condition Eq. <I5> requires that the electric field 
in the comoving frame at the inner radius vanishes: v,, || Bp — » Zs^ = 
and Ep = ClRi„Bp/a(ri„). The condition on implies that the ra- 
dial component of the magnetic field remains constant. We chose 
for B,(ri„) different values to investigate both cases with low mag- 
netizations and high magnetization. Bg and B^ where extrapolated 
using 0-th order reconstruction (we found that linear interpolation 
can lead to spurious oscillations). The value of the tangential ve- 
locities were derived using the remaining constraint of the frozen-in 
condition: vg = v^Bg/Br and = SlRj„/a(rj„)-v,.B^/Br. Given that 
the 3 components of the velocity are derived independently there is 
no guarantee that v' < I, so care has to be taken to avoid sharp 
transients and spurious oscillations in the tangential magnetic field 
near the inner boundary. 

In Table the parameters of our various models are given. 
Equations <17> - <22> show that the problem can be parameterized 
in terms of the ratio <b- IT (assuming B^ scales as Br) and not on 



the specific value of density and magnetic field; more generally the 
parameter governing the properties of the system is ctq = (^^Qr/T 
{ Michel 1969). The bulk of our simulations have been done using a 
fixed value for Q., in order to allow a more straightforward compar- 
ison among the various results; however in a few cases (BI of the 
2D monopole, and B I , B2 of the 2D dipole) we use a different rota- 
tion rate to check whether the energy and angular momentum loss 
rates indeed only scale with ctq. We note that it is computationally 
more efficient to increase the magnetic flux than to drop the mass 
flux in order to achieve higher ctq. Only in the ID case, where reso- 
lution is not a constraint, are we able to investigate the behavior of 
the system for different values of <S. 



5 RESULTS 

5.1 The ID Monopole 

As a starting point for our investigation, we consider the simple 
case of a relativistic monopolar wind in ID. This is the relativis- 
tic ex tension of the classic W eber-Davis solution for a magnetized 
wind jWeber&Davisll967h . and represents a simplified model for 
the flow in the equatorial plane. The ID model can also be used 
both to verify the accuracy of the code and as a guide in understand- 
ing the 2D simulations ( ^5.211531 . Here we assume vg = Bg = 0. 
The solenoidal condition on the poloidal magnetic field reduces to 
B, oc r^. 

We calculate a number of wind solutions for different values 
of the flux O and the entropy <S. The results are summarized in 
Table |2| Case AO is a reference case for an almost unmagnetized 
wind, and it will be used for comparison with the weakly magne- 
tized regime (A~D). As stated above, the solution can be parame- 
terized in terms of <to = <t>-Q?/T'. We have verified that the mass 
loss rate T depends strongly on the value of the sound speed at r,„, 
and drops rapidly as the pressure approaches the critical value for 
H > 0, which also depends on the magnetic field strength. In con- 
trast, over the range of parameters studied, f has a relatively minor 
dependence on the value of the magnetic field at fixed <S: increasing 
Br(ri„y by 3 orders of magnitude corresponds to an increase in f 
by just a factor ~ 1.7 (compare models AO and FinTab.|5J. The rea- 
son for this is that in all cases listed in Table|5|the AL point is larger 
than the SM point. Thus, the magneto-centrifugal effect of increas- 
ing the density scale height in the region interior to the SM point, 
where the mass loss rate is set, is already maximized. For all the 
cases investigated the value of Tp/p(Q,r)^ at the AL point is always 
less than 0. 1 (in case A it is ^ 0. 1 while in case E it is ^ 0.02). All 
the solution can thus be considered centrifugally driven. We expect 
and find a sharp drop in !F as we go to yet smaller magnetization. 
For example, in the purely hydrodynamical case without magnetic 
fields, T a 30, six times lower than !F for Case F. 

Note that the value of y reported in Table|5|is at a fixed po- 
sition, r = lOOrf^s, which is generally outside the FM point (see 
Fig.|2}. Case I is an exception. For reference, the Lorentz factor at 
~ 300 rjfs is 7.8 for this model. 

All values in Table |2| are given in code units. The following 
relations can be used to scale to physical units, in terms of the value 
of density and magnetic field at the injection radius r,„. The mass 
loss rate is 
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Table 1. Parameters of the numerical models (unit with c = 1). 



1 T~\ TV fl^ J 

ID Mod. 


AO 


A 


B 


C 


D 


E 


F G H 1 


t^vK^in) !p\Xm) 


0.033 

U.UU4 


0.033 


0.028 

0.0*4 


0.025 


0.021 

0.0*4 


0.018 
o.o^ 


0.033 0.030 0.030 0.030 

4 o jZ oO 


2D Monop. Mod. 


A 


B 


Bl 


c 


D 


E 




Br(rin9lp(rm) 


0.033 
0.04 


0.033 
0.4 


0.033 
0.4 


0.033 
4 


0.033 
40 


0.033 
200 




2D Dip. Mod. 


A 


B 


Bl 


B2 


c 






BAr,„,e = oflp(n„) 


0.033 
0.64 


0.033 
6.4 


0.033 
6.4 


0.033 
6.4 


0.033 
64 







In all models £1 = 0.143 except models Bl which have SI = 0.214 and model B2 which has £1 = 0.0715. See Eq. 1241 for the corresponding val ue of the 
rotation period. c\. = Ypjiph). Case AO is a reference case for an almost thermally driven wind. In term of the standard wind parametrization ^S akurai 
^985, Daiene & Drenkhahn 2002) all our cases are centrifugally driven: the value ofTplp{Q.r)^ at the AL surfaces, is always less than 0.1 except in case 
AO where it is 0.5. The value of the lapse at the injection radius is a(rj„) = 0.79 coiTespondig to an escape speed 0.6c. The unit of length corresponds to 
the radius of the neutron star r^s , the unit of time is r^s /c 



the rotation period is 

Ulkm/ 

the surface magnetic field strength is 

M 



(24) 



B.(r,„) = 4.25- 10' 
P 



" ' 1.9 ■ W-'*MeS-' 
P(r,n) ' - 



\1.6-10-3s/ llOiOgcm-V Ulkm/ 



^ID'Ogcm- 
the total energy loss rate is 



(25) 



6.95 • 10" 



B,(R,„) 
\(D2n2A8.7- 10" G 

P 



Ulkm/ U.6 • 10-3s/ 



.6 • 10-3s 
the angular momentum loss rate is 



ergs s 



(26) 



1.57 ■ 10* 



\(D2£lll 8.7 • 1013 G 



\llkm/ 



(l.6- 10-3s) 



ergs, 



(27) 



where J = £T and H = 'HT. 

In Fig.|2|we plot the velocity profiles for cases F, G, H, and 
I, together with the location of the SM, AL, and FM points. Note 
that the position of the SM point does not change significantly and 
is given roughly by Eq. ^5). In Fig.|3|the angular momentum loss 
rate and energy loss rate are plotted as a function of ctq. The con- 
vergence to the force-free solution is evident (see also Table|2}. An 
alternative way to parametrize how close the soution is to the force- 
free limit is by considering fa/Rl iDaigne & Drenkhahn 2002*), as 
shown in Tab|5| In the mass-loaded cases (ctq < 1; models A-D) 
we find that J/(t>^Cl oc M'/^ in accord with the expectations from 
82.21 We find that all the points (we excluded case AO because we 
are interested in case where > rsM, typical of the transition from 
mass loaded to force-free) can be approximated by a relation of the 
form j l<SrQ. = Co + Ci(l/cro)''^, where Co should correspond to 
the force-free limit J/(t>^Sl = 1. A fit to our results gives Co = 0.98 
and Ci = 0.93. A similar expression can be written for the energy 
losses. At low (To the total energy loss rate scales as the mass loss 



Table 2. Results of the ID models. 


Mod. 


o"o 


T 






yilOOi-Ns) 


rAlRt 


AO 


0.013 


101 


10.6 


99.0 


1.12 


0.30 


A 


0.084 


142 


4.17 


18.0 


1.20 


0.46 


B 


0.164 


74 


3.34 


9.98 


1.22 


0.56 


C 


0.307 


39 


2.71 


6.06 


1.32 


0.65 


D 


0.634 


19 


2.18 


3.75 


1.42 


0.75 


E 


2.81 


4.25 


1.52 


1.86 


1.87 


0.89 


F 


6.85 


175 


1.31 


1.49 


2.30 


0.93 


G 


18.4 


130 


1.17 


1.24 


2.94 


0.96 


H 


61.4 


155 


1.08 


1.10 


4.00 


0.98 


I 


127 


188 


1.06 


1.06 


4.84 


0.99 



j = JiT and H = l-iT . Values of T are given in code units, see Eqs. 1231 - 
1271 for conversion to physical units. In all cases n = 0.143. 



rate, as expected, if/O^fi^ oc M. For large ctq, the solution con- 
verges to the force-free limit, Hl<b-Q} = 1 (see Fig.j^. We find 
that the transition between the two limits can be fit with H/9^n = 
Co + Ci( 1/0-0)°'', where C„ = 0.98 and C, =2.16. We also con- 
sider the reduced energy loss rate (the difference between the total 
energy and the rest mass energy), and find that this too can be ap- 
proximated with a power law: (H - T)/<^^Q? = Co + Ci(l/o-o)"^', 
with Co = 0.98 and Ci = 0.965. We stress that these power law 
relations have been determined by fitting the results of our sim- 
ulations. As such, there is no guarantee that these trends can be 
extended outside the range we investigated, but they do correspond 
well to the estimates given in ^ 

Our simulations focus on the region close to the neutron star 
and so the problem of the acceleration of the outflow at large dis- 
tances ca nnot be properly addressed. However, as shown, for ex- 
ample, bv lPaigne & Drenkh ahn ( 2002) the efficiency of conversion 
from magnetic to kinetic energy in the strict monopole limit is very 
low. Faster than radial divergence in the flux tubes is required af- 
ter the FM point to increase the acceleration significantly ( 82. U . 
Within our computational domain in Case I we find that the ratio 
of particle kinetic energy flux to electromagnetic energy flux scales 
approximatively as r''^ and shows a tendency towards saturation. 
At a radius of 300r;vs in Case I, the ratio is still less than 5%. 
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Radial velocity 



Table 3. Results of the 2D monopole models. 




100 



1000 



(km) 



Figure 2. Radial velocity, and position of the SM (plus), AL (diamond) and 
FM (triangle) points, in the ID monopole case. From bottom to top lines 
refer to cases F. G, H and I. 



Loss rates 




0.001 0.010 0.100 I.OOO 10,000 100.000 

Figure 3. Loss rates for the ID monopole calculations in non-dimensional 
units (Table l2l. Upper panel: angular' momentum loss rate. Lower panel: 
total energy loss rate. Dashed curves represent the theoretical expectation 
for the losses in the mass loaded cases j oc M^l^ and H oc M. Continuous 
curves represent the best power-low fit given in the text. Dotted lines are the 
force-free solution. 



5.2 The 2D Monopole 

The ID monopole discussed above does not take into account de- 
formations of the poloidal field lines by the moving plasma. As a 
consequence, the conversion of magnetic energy into kinetic en- 
ergy of the accelerated wind is inefficient. To understand if and 
how deviations from a strict monopole may affect the dynamics of 
the outflow it is necessary to perform 2D simulations. We focus 
only on the region close to the star, within 100 - 200 rt^s , and we 
consider both mass-loaded (cto < 1) and Poynting-flux dominated 
((To > 1) regimes (see Tab les[T]&|3. In contrast to the cases consid- 
ered by Boeovalov (2001), where the mass flux was fixed at r„, and 
a cold wind (p = 0) was assumed, here the mass flux is derived self- 
consistently, with pressure at the base of the wind being the control 
parameter for the flow. Even if Cjir,,,) ~ 0.1c, the dilference with 
the pressureless case is not trivial. For example the location of the 



Mod. 


D"0 


T, 








A 


0.121 


98 


1.23 


11.9 


0.40 


B 


1.00 


119 


1.11 


2.28 


0.68 


Bl 


1.22 


217 


1.05 


1.99 


0.70 


C 


9.67 


122 


0.755 


0.875 


0.92 


D 


68.5 


173 


0.678 


0.699 


0.98 


E 


211.5 


280 


0671 


0.679 


0.99 



{4jT)-'_[rds and (Tq = D^a^/r,. Values are given in code units, see 
Eqs. I23I - I27I for conversion to physical units. The value of .S in all cases is 
0.018. In all cases il = 0.143, except case Bl which has Q. = 0.214. rAgguat 
is the radial distance of the AL sulface on the equatorial plane. 

FM point in the ID monopole is at infinity if p = 0, so, in princi- 
ple, one might expect a higher efficiency also in the 2D case. More 
important, in our case, the velocity at the base is much smaller that 
c, so that collimation in the region close to the star could be more 
efficient. 



5.2.1 Magnetic, Mass-Loaded Winds 

In Fig.|4]we show the results for a heavily mass-loaded case (Case 
A) corresponding to a ctq = (47r)"' ^, IT ds = 0.121, where 
the integration is performed over 4n solid angle. 

At Tin the mass flux profile scales approximatively as sin^ 9, 
and is minimal at the pole. As a result of magnetic acceleration and 
centrifugal support at the equator, the mass flux is higher than in the 
corresponding ID hydrodynamical non-rotating case and it is about 
the same as in the ID monopole. However, at the pole the mass 
flux is lower because of magnetic collimation on the axis (Kopp & 
Holzer 1976). The large difference in mass flux between pole and 
equator is manifest in the elongated shape of the FM surface. The 
upper left panel of Fig.|4|shows that the field lines are very close to 
radial and that scales approximatively as sin 6. At the pole, the 
FM surface falls outside of the computational domain, whereas the 
FM surface intersects the equatorial plane at S.lr^s (Xin = 1-1''a's)- 

The SM surface also has a large axis ratio: it intersects the 
pole at a distance of llrjvs and it intersects the equator at LSr^/^ 
It is interesting to look also at the position of the AL surface. Its 
distance from the Light Cylinder is an indicator of how close the 
solution is to the force-free limit and it strongly reflects the degree 
of magnetization. While on the axis where = the AL and 
FM surfaces are coincident, away from the pole the toroidal field 
component does not vanish and the two critical curves separate. 
The AL surface intersects the equator at a distance of about 2.1 r^s 
(to be compared with Rl = G.Sr^s )• 

One might expect the Lorentz factor to be largest in the equa- 
torial region, as a result of stronger magnetocentrifugal effects 
there. However, contrary to this expectation, we observe that y 
peaks at about 70°. Such a result was also obtained by Boeovaloy 
(2001) for cold flows. This effect is stronger in our calculations 
because of the lower overall Lorentz factor. We also notice the ex- 
istence of a very slow channel along the axis. We want to stress 
that the FM surface is outside the computational domain within 3° 
of the axis, and so the solution has not converged fully in this re- 
gion. However, by increasing the radial computational domain, we 
find that the main effect of failing to capture the FM surface at the 
pole is that the wind in this region is less coUimated and somewhat 
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Figure 4. Results for the 2D monopole in the weakly magnetized Case A (Tablel3l. Upper left: contours represent poloidal magnetic field lines, while colors 
represent the ratio IB^/B, !. Upper right: colors and contour represent the Lorentz factor. Notice the presence of a slow channel on the axis and the peak in 
velocity at about 70°. Lower left: angular momentum flux in adimensiona units £.T H'\n<S)^Q.). Lower right: total energy flux fiT KAn^^Q?') in adimensional 
units (see Eqs. I23I - I27I for conversion in physical units). Note that these fluxes peak at high latitudes (compare with Fig.lsl. 



faster than it should be. So we expect the wind to be more col- 
limated and slower, with yet larger computational grids. Whether 
the FM surface is inside the computational domain at the pole or 
not has relatively little effect on those streamlines at lower latitudes 
that do pass through the FM point. Typical deviations are found to 
be less than 1%. 

As the lower panels of Fig. |4| show, we find that both the en- 
ergy flux (which is mainly kinetic) and the angular momentum flux 
peak at high latitudes. In addition, the mass flux at large distance 
from the neutron star is higher close to the axis because of magnetic 



coUimation. In contrast, the conversion of electromagnetic energy 
to kinetic energy is maximal along the equator, even though in this 
mass-loaded Iow-ctq case the electromagnetic energy flux is lower 
than kinetic energy flux, and the wind terminal Lorentz factor is 
mainly given by the conversion of internal energy to kinetic energy. 
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Figure 5. Results for the 2D monopole in the highly magnetized Case E (Table|3). Upper left: contours represent poloidal magnetic field lines, while colors 
represent the ratio \B^/Br\. The FM surface (dotted line) is more distant from the axis while the AL surface (dashed line) is very close to the Light Cylinder. 
Upper right: colors and contours represent the Lorentz factor. There is no evidence here of the relatively slow channel on the axis as in Case A (see Fig.|3 and 
the Lorentz factor scales as sin (9). Lower left: angular momentum flux in adimensiona units X?^/(4ff<I'^n). Lower right: total energy flux'Kr/(4ffa)2n2) in 
adimensional units (see Eqs. I2ji - I27l for conversion in physical units). Here, these fluxes are higher at the equator than at the pole (compare with Fig.|3, as 
expected when the flow is relativistic and Poynting-flux dominated. 



5.2.2 Poynting-Flux Dominated Winds 

In Fig. 13 we show the results for a Poynting-flux dominated flow 
(Case E), corresponding to a ratio ctq = (4;r)"' (b^Q?- IT ds = 
21 L5, where the integration is again performed over 4n sohd angle. 
In this case the solution is mostly magnetically driven, and plasma 
effects lead to small deviations with respect to the force-free limit. 

Similar to the mass-loaded Case A, at the neutron star sur- 



face we find that the mass flux is higher at the equator than at the 
pole, and that the AL and FM surfaces are extended in the direction 
of the pole. Here, the FM and AL surfaces intersect the equatorial 
plane at a distance of AGr^s (to be compared with 40.5 = cTq^RiJ 
and 6.1rf^s, respectively, while the SM surface is more spherical 
than Case A and intersects the pole and the equator at a distance 
of 2.5rNs and Llr^s, respectively. Magnetic field lines are again 
monopolar, but while in the mass-loaded Case A this was a conse- 
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quence of a amrginally centrifugally driven wind, now this is due 
to electromagnetic force balance. The dynamics of a magnetized 
outflow are governed by the combination of Lorentz and Coulomb 
forces. In the relativistic regime the Coulomb force cannot be ne- 
glected and, as the flow speed approaches c, the Coulomb force 
balances the Lorentz force, suppressing collimation. As a conse- 
quence, the flux tubes in a Poynting-flux dominated wind have an 
areal cross-section that scales as r-. It is known that the efficiency of 
conversion of magnetic energy to kinetic energy incr eases when the 
flow divergence becomes more than radial (Daigne & DrenkhahnI 
l2002h . In the high-cro simulation presented here, at least within the 
limited computational domain employed, we do not s ee evidence 
for ef ficient conversion and acceleration. As found by iBoeovalo^ 
l200lh there is evidence for a narrow coUimated channel very close 
to the axis, but in our case, where the mass flux at the injection is 
not imposed, the mass flux in the channel is not strongly enhanced. 
However the FM surface on the axis does not close in our compu- 
tational box so no strong conclusion can be drawn. 

The upper right panel of Fig.|5|shows that the Lorentz factor 
in the wind scales approximatively as sin0. The maximum value 
achieved in the computational domain is 8. The latitude dependence 
is not appreciable, probably because our solution does not extend 
far enough away from the star. As in the ID models, for the range 
of parameters investigated, the total mass flux from the star is not 
much affected by the value of the surface magnetic field. For the 
Poynting-flux dominated Case E it is about 3 times higher than in 
the mass-loaded Case A, despite the fact that the magnetic energy 
density is 5000 times higher. This again follows from the fact that 
the AL surface is larger than the SM surface. However, there are 
important qualitative differences between Case A and Case E. In 
case E we find that on large scales magnetic acceleration is domi- 
nant and the mass flux is maximal on the equator, instead of on the 
axis. In addition, the lower panels of Fig.|5|show that the energy and 
angular momentum fluxes increase toward the equatorial plane and 
are almost completely magnetically-dominated. In fact, the energy 
and angular momentum loss rates scale as sin^ 6, as expected in the 
force-free limit. This strong transition in the angular dependence 
of the energy and momentum fluxes, from strongly collimated to 
equatorial, occurs at (Tq ~ 1. 

Figure |5| shows the behavior of the angular momentum and 
energy losses as ctq increases. We again observe the convergence 
to the force-free limit (see also Tab. |3} and recover the expected 
behavior J/<^^C1 oc M''^ in the mass-loaded cases. The conver- 
gence to the force-free solution can be approximated, as in the ID 
case, with a power law of the form Co + Ci{llcrof. The fits are 
different than the ID models because the flux tubes deviate from 
being strictly monopolar. We find for the angular momentum loss 
rate 7/<l)^n = 0.664 -l- 0.42(l/a-o)°'^', for the total energy losses 
///O^n^ = 0.660 + 1.66(l/o-o)''•*^ and for the reduced energy 
(H - r)l<i>-il- = 0.66 + 0.60(l/o-o)°-^^ As mentioned before, we 
have derived these fits within the parameter range of our simula- 
tions, and there is no guarantee that they can be extrapolated for 
much higher or lower magnetizations. We stress that in the highly- 
magnetized Case E our results are v ery clos e to the force-free so- 
lution: y/O-fl = H/<t>^Q.- = 2/3 (Michel 1991). Note also that 
such a value is lower than what would be expected from a trivial 
2D extension of the ID monopole (OJS=J^'' sin (df d6). 

The efficiency of conversion of the magnetic energy to kinetic 
energy is maximal on the equator, but it does not exceed 10% at 
the outer boundary. We note that conversion is faster than logarith- 
mic in radius and at the edge of the computational box it scales 
as similar to our ID results. However, we cannot draw strong 
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Figure 6. Loss rates for a 2D monopole in non-dimensional units. Upper 
panel: angular momentum loss rate. Lower panel: total energy loss rate. 
Dashed curves represent the theoretical expectation for the losses in the 
mass loaded cases J oc M'/3 and H oc M. Continuous curves rep- 

resent the best power-low fit given in the text. The dotted lines are the 
force-free limits. The squai'e mark indicates Case Bl, which has a differ- 
ent rotation rate (Tab.fT]. 

conclusions on the terminal efficiency far outside of the FM surface 
because of the limited size of our computational domain. As shown 
by Bogovalov (2001), y seems to increase after the FM point and 
then saturates at larger distances. 

5.3 The Aligned Dipole 

To lowest order, currents in the neutron star should generate a 
dipole magnetic field. This field configuration is much more re- 
alistic than the monopolar models considered in ^5. H and ^5.21 A 
dipole field may also have interesting consequences for the asymp- 
totic character of the outflow. For example, it is possible that the 
presence of a closed zone, outside of which the open field lines at 
first expand much more rapidly than radially, might provide for a 
more efficient conversion of magnetic energy into kinetic energy, 
leading to a higher terminal Lorentz factor (0. Below, we review 
several numerical issues associated with our dipole wind solutions 
and then we present our results, which again bridge the transition 
from low- to high-cro outflows. 

5.3.1 Numerical Challenges 

The modeling of a magnetic wind with a closed zone and an equa- 
torial current sheet presents a number of numerical difficulties. We 
have encountered two problems in particular that bear mention. The 
first is that the outer edge of the closed zone rotates/a.sfer than what 
is required by Eq. J19> . We believe this "super-corotation" is con- 
nected with numerical dissipation in the equation for the evolution 
of and that it stems from the fact that the boundary between the 
closed and open field regions is not grid-aligned.'* By suppressing 
the upwinding term in the HLL flux we were able to reduce the 
deviation from the corotation condition from 10 - 20% to ~5%, 
but at the price of making the code less stable when the flow is 

* hicreasing the resolution of the simulation or using a chai'acteristics- 
based solver iKomissaro\l 1 99^1 private communication) does not improve 
the accuracy of the solution. 
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Table 4. Results of the dipole models. 



Mod. o-o r, J/lija H/O^Qr ry/Rt 



A 


0.298 


66 


1.53 


5.28 


0.26 


B2 


1.34 


7.8 


0.965 


1.81 


0.32 


B 


2.77 


43 


0.853 


1.27 


0.37 


Bl 


3.91 


114 


0.816 


1.10 


0.48 


C 


17.5 


42 


0.701 


0.73 


0.47 



(4g)~ ' Tds and ctq = OjSl^/'Fi. Values are given in code units, see 
Eqs. <231 - <27t for conversion to physical units. Cases Bl and B2 have dif- 
ferent rotation rates (Tab. IT]. All cases have S = 0.018. 

highly magnetized. Unfortunately, of the previous papers dealing 
with winds in the presence of a cl osed zone in the MHD regime, 
onlv lKeppens & Goe dbloecf l'2000') discuss deviations from corota- 
tion in the closed zone. In their paper these amount to ~ 10-20%, 
and they consider only parameters appropriate to the Sun, a slow 
rotator. 

The second problem is that the magnetic field undergoes re- 
connection at the neutral current sheet on the equator close to 
the position of the Y point, where the last closed field line inter- 
sects the equatorial plane. As a consequence, plasmoids are formed 
and advected away, thus preventing the system from reaching a 
steady-state configuration.^ At the current sheet and change 
sign, the MHD approximation fails, and a sharp discontinuity de- 
velops that cannot be well-resolved. Figure |7| shows an example 
of plasmoid formation. Although it is well known that current 
sheets in the presence of a K-type point are subject to reconnection 
and the continuous formation of plasmo ids j&ideve & Leer 20031 
iTanuma & Shib ata 2005 ." Yin et alj2000l) . in our case the high value 
of (To does not allow us to properly resolve the current sheet. For 
this reason, the reconnection processes are dominated by the intrin- 
sic numerical resistivity of our numerical scheme. 

We have found that the formation and growth of plasmoids 
depends on t wo terms in the d efinition of the HLL electric field 
(see Eq. 44 of iDel Zanna et al.li2003.') : dBr/de and B,-ve. The first 
behaves like an explicit resistivity and seems responsible for the 
evolution of the plasmoids as they are advected off the grid along 
the equator. The second term behaves like forced reconnection and 
controls the initiation of plasmoid formation (given that the current 
sheet is not resolved, vg is not reconstructed to zero on the equa- 
tor). Setting both terms to zero causes the closed zone to disappear 
entirely and the system evolves toward a modified split monopole. 
To deal with this issue and explicitly enforce a steady state, we 
opt for the following procedure: in a calculation with plasmoids we 
note the position of the Y point and we then impose the conditions 
dB, /d6 = and B^vg = on the equator outside the position of 
the Y point as inferred from the calculation without these bound- 
ary conditions. The condition dBJdO = can be justified because 
of the structure of the characteristic waves and because the solu- 
tion should be symmetric about the equatorial plane. The equator 
is a contact discontinuity, so only the total pressure is important 
as a boundary conditions and not the sign of the parallel magnetic 
field, in fact for infinite conductivity a pure monopole and a split 



^ Wind calculations performed over a full 180 degrees show that the for- 
mation of plasmoids is not an artifact of our 90° computational domain. 
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Figure 7. Magnetic structure of a magnetically dominated flow (Case C). 
The upper panel shows the poloidal magnetic field structure, with a snapshot 
of outflowing plasmoids fonning along the equatorial current sheet outside 
the closed zone. The lower panel shows a blow up of these plasmoids, which 
travel out at close to the speed of light. The color indicates the angular 
momentum density, with roughly 20% of the angular momentum loss being 
canied by these intermittent structures. Values are in code units (Eqs. 1231 - 
1271 '). As stated in the text, the numerical values in the plasmoids depend on 
numerical resistivity. 
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monopole have the same solution. The condition B,vg = corre- 
sponds to enforcing Vg = 0. 

We note that the steady state solution need not be the physi- 
cally correct solution found in nature. The equatorial current sheet 
in the vicinity of a neutron star is undoubtedly dissipative, and is 
thus subject to reconnection and plasmoid formation. However, a 
full understanding of this behavior requires — at the very least — 
the use of resistive RMHD, so that the dissipation can be controlled, 
rather than being fully numerical, as in our current calculations. 
Such a treatment is beyond the scope of this paper. Thus, we focus 
our attention here on the forced steady state solutions described 
above. As a test, we have compared the global energy and angu- 
lar momentum loss rates between time-dependent calculations with 
plasmoids and those with our forced steady state boundary condi- 
tions outside the Y point. In general, losses are higher in the latter 
calculations because the open magnetic flux is larger. In Case A the 
difference is less than 5%, while in Case C it is about 15-20% (see 
Tab.|31and Fig.[S). In the time-dependent calculations the individ- 
ual plasmoids represent fractional deviations in j and H from the 
average of up to 15% in Case C, and less for lower ctq flows. Thus, 
the plasmoids do not appear to be that dynamically significant for 
the overall energy and angular momentum losses from the neutron 
star. 
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Figure 8. Loss rates for a 2D dipole in non-dimensional units. Upper panel: 
angular momentum loss rate. Lower panel: total energy loss rate. Dashed 
curves represent the theoretical expectation for the losses in the mass loaded 
cases j IX M''^ and H oc M. Continuous curves represent the best power- 
law fit of the 2D monopole of Fig.|6] The dotted lines are the force-free 
solution. The squares mark cases Bl and B2, which have difl^'erent rotation 
rates (see Tab.fH. 



5.3.2 Results 

In Table |4| we show the results of our simulations using the pre- 
scription described above. In the monopole models the magnetic 
field lines are open so that normalizing to the surface magnetic field 
or to the total magnetic flux O is equivalent. In the dipole case, be- 
cause of the closed zone, the two quantities are not proportional. 
We find it useful to normalize in terms of the magnetic flux eval- 
uated on open field lines, Oo, where equations jl7> - J22t hold. As 
a consequence, we can define an equivalent surface magnetic field 
Br-eqmv(ri„) s& the surface magnetic field of a monopole that has 
the same amount of open magnetic flux. Contrary to the monopole 
case, increasing the magnetic field strength at the stellar surface 
reduces the mass loss rate because the size of the closed zone in- 
creases. Similarly, for the parameters explored here, an increase in 
the magnetic field strength by a factor of, say, ~ 3 at the stellar sur- 
face leads to an increase in the open magnetic flux 0„ by a factor 
of just ~ 2.5, rather than the one-to-one behavior for the monopole. 
We can derive an approximate relation between the magnetic field 
at the pole on the surface of the star {r^s ), and the equivalent sur- 
face magnetic field, that can easily be computed if H and Q. are 
known. We find that in all our cases the relation is 



Br(rNs,e = 0) ^ 1.6S,,_,„,„v(rN5) X (rylrNs)- 



(28) 



Thus, the problem of relating the surface magnetic field to the loss 
rates is reduced to the problem of determining the size of the closed 
zone. 

In Fig.|8|the angular momentum and energy loss rates are plot- 
ted. The continuous lines in this figure are not a fit to our dipole 
results, but are simply the same curves as in Fig. |6|/or the 2D 
monopole. In addition, as in Fig.|S| the dashed lines are the ana- 
lytic expectation for the loss rates in the mass-loaded limit /or the 
monopole ( 32. 2t . This shows that if the solutions are parameterized 
in terms of the open magnetic flux, then the dipole and monopole 
winds have very similar behavior within the parameter space we 
have investigated. This can be easily understood if one considers 
the structure of the outflow in the far region. Given that H and / are 
integrals, they can be evaluated at any distance from the star. Even if 



the the field is dipolar close to the star it is nearly monopolar outside 
Rf Even when we do not impose our steady-state boundary condi- 
tions at the equator, and plasmoids are present during the evolution, 
H and j closely follow our results for the monopole (again, as long 
as these losses are written in terms of the open magnetic flux). Note 
also that in Figure|8| we find that H/^lCl- and J/^lSl con verge to 
the ex pected value of 2/3 in the force-free limit I Conto poulos et alJ 
19991 lGruzinovll200.5h . However, as we discuss below, ry < Rl in 
all of our calculations, contrary to the assumption that ry = Rl in 
the above force-free treatments. 

We can extrapolate our results to the force-free limit for the 
spin-down rate. In terms of the equivalent surface magnetic field, 
and using Eq. <28> with ry = Rl and 6,, = Q.5B, {rNs , 6* = 0) we 
have: 

2 

NS ^r—equiv 



H 



24 , 4 
- 25^^" 



(29) 



in agreement with lGruzino\i iioO^ . 

The flow structure for both low ctq (Case A) and high ctq (Case 
C) are shown in Figure|5| These figures show that the poloidal mag- 
netic field for R > Rl has a structure very similar to that obtained 
from the monopole calculations. Also, B^ scales as sin 0, except in 
a region close to the equatorial plane where it changes sign. This 
agrees with the results of Fig.|8|showing that the energy and angular 
momentum losses scale as for a monopole. 

As expected, we find that the size of the closed region in- 
creases with magnetization. The position of the Y point moves from 
l.Srivs in Case A, to 2.5r^s in Case B, to ^.Ir^s in Case C. Thus, 
an increase in (Tq of a factor of ~60 corresponds to a ~70% in- 
crease in the size of the closed zone. Importantly, even at relatively 
high (To (~ 18 for Case C), these values are significantly smaller 
than the Light Cylinder radius Rl = 6.8rjv5 . In order to understand 
the systematics of the Y point, we have calculated several models 
with different rotation rates CI (see Tabs.Q&|4}. The radius of the 
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Figure 9. Result of the 2D dipole. Upper panels represent case A, lower panels case C . Left: contours represents poloidal magnetic field lines, while colors 
represent the ratio B^/B,-. The dotted line is the FM surface, the dashed line the AL surface. In case C the AL surface is close to R^. Middle: contours represents 
poloidal magnetic field lines, while colors represent the ratio B^/B^, in a region closer to the star. Now the dotted line represents the SM surface. Right: colors 
and contour represent the Lorentz factor. 



Y point changes from l.lr^s in Case Bl, to l.Sr^s in Case B, to 
4.3rf^s in Case B2. Although ry is largest in the model with slowest 
rotation (Case B2), this model has the smallest ratio ry/Ri^. Indeed, 
we find that ry/Rt decreases as CI decreases. This trend yields a 
braking index less than 3. If one takes Cases Bl, B, and B2 as a 
time-series in the life of a neutron star, one would infer a braking 
index ~2.2. 

In general the size of the closed zone will depend on the phys- 
ical conditions at the stellar surface, including the thermal sound 
speed and the mass density, which govern the mass loss rate on 
each open streamline (see jmiMestel & SpruiiigSTh . This shows 
that even a small thermal pressure can have an important effect on 
the torque. Figure|9|shows that in Case A both the AL and the FM 
surfaces are inside the radius of the Light Cylinder. However, in the 
high-cTQ Case C the AL surface is very close to R^. This suggests 
that although the AL surface rapidly approaches Rt as ctq increases, 
the position of the Y point remains inside Ri and is a weak func- 
tion of (To. It is thus possible to produce a relativistic outflow, even 
if the Y point is well inside Ri^. Because the range of parameters 
we have explored is fairly limited, we can only conclude that quite 
large magnetization is required to achieve ry = R^. A rough extrap- 
olation of our results implies that B,(rj„,6 = 0)^/p(r;„) must be of 
order 10^ - 10* to achieve ry « R^, but at the moment we deem 
premature to draw strong conclusions. 

We stress here that the magnetization parameter ctq is defined 



as an integral average. It is known that in a dipolar field geometry 
with sub-slow magnetosonic injection, the mass flux at the edge of 
the closed zone is higher than the average integrated value over the 
entire star (Kodd & Holzet 1976). Thus, the magnetization varies 
from high to low latitudes. For example, in our Case C the magne- 
tization on the open flux tubes immediately adjacent to the closed 
zone is less than 7, while the ctq parameter for this model is 17.5. 
Because the position of the Y point and the shape of the magneto- 
sphere depend on local equilibrium between the closed and open 
zones, this strong variation in magnetization in latitude may help 
explain why even if the global flow has ctq » I, ry is less than Rt. 

The upper and lower far-right panels of Fig. |5| show that for 
the Lorentz factor we recover similar behavior as for the monopole 
(compare with Figs.|31&|5j. In Case A, the Lorentz factor reaches 
its maximum at about 70° from the equator. We also notice that 
there is now a slower equatorial flow corresponding to what is 
known as the "slow solar wind" in models of the Sun's outflow. In 
Case C the Lorentz factor appears to scale mostly as the cylindri- 
cal radius. This again was found for the monopole. The slow wind 
region is still present, but now the boundary conditions we have 
imposed on the equator to suppress the plasmoids cause a noisy 
structure in v,-. 

The shape of the field lines outside the FM surface does show 
significant diff'erences between the monopole and the dipole. In the 
dipolar Case C at r = 50rjv5 the maximum value of Bg/Br = -0.06. 
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Such a value for this ratio is in between Cases C (-0.08) and D 
(-0.03) for the 2D monopole. Our limited computational domain 
prohibits a more quantitative study of the poloidal field line shape 
at still larger distances from the neutron star. 

Lastly, the angular distribution of the energy and momentum 
loss rates in our dipole models are qualitatively similar to those we 
obtained for the 2D monopole (see the bottom panels of Figs.|31 
&|5)- For low-cTo outflows, because of hoop stress, H and / far 
from the neutron star are peaked at high latitudes, along the axis of 
rotation. In contrast, for high-cTo outflows, H and J are maximum 
near the equator, similar to Fig.|5|(Case E). 



magnetization along the open field lines nearest the closed zone is 
less than 7. 

(vii) Over limited dynamic range, the ratio rylRt decreases as 
0) decreases in the dipole models. This behavior is consistent with 
a braking index less than 3 (see ^2. H and ^6.H . 

Concerning the acceleration of the wind, we do not find any evi- 
dence for efficient magnetic to kinetic energy conversion. However 
we want to stress again that our computational box does not ex- 
tend far enough from the star to appreciably probe the asymptotic 
behavior of the wind. 



6 DISCUSSION & CONCLUSIONS 

We have solved for the dynamics of time-dependent relativistic 
MHD winds from rotating neutron stars, including the efl'ects of 
general relativity on the gravitational force. The mass loss rate 
along open field lines is derived self-consistently as a solution to 
the RMHD equations, subject to the boundary conditions at the 
stellar surface (finite thermal pressure) and a polytropic equation 
of state. We consider ID and 2D monopole field geometries and 
the aligned dipole, and we explore solutions that cover the param- 
eter regime from non-relativistic mass-loaded Iow-ctq outflows to 
relativistic Poynting flux dominated high-cro outflows. Our primary 
results are: 

(i) In the ID and 2D monopole calculations, we reproduce the 
expected analytic trends in both the high-cro and Iow-ctq limits. 
In particular, the solutions asymptotically approach the force-free 
limit when cto s> I . 

(ii) In both the dipole and 2D monopole solutions, when ctq < 
1, the energy and momentum losses far from the neutron star are 
highly directed along the axis of rotation. The zenith angle at which 
these fluxes are maximized is an increasing function of ctq so that 
for (To > 1, the losses are primarily equatorial. 

(iii) In both the dipole and 2D monopole solutions, for winds 
with (To < I, the Lorentz factor peaks at high latitudes. 

(iv) For the aligned dipole, the equatorial current sheet may be 
unstable to the formation of plasmoids, leading to time-dependent 
spindown of the neutron star. A proper treatment of dissipative 
RMHD in the equatorial region is needed to explore this issue more 
completely. 

(v) If the energy and angular momentum losses from the aligned 
dipole are parameterized in terms of the open magnetic flux, then 
the results are nearly identical to those from the 2D monopole 
solutions in both the mass-loaded and Poynting flux-dominated 
regimes. In particular, the dipole calculations quantitatively ap- 
proach the force-free limit when (To » 1. The normalization for 
the angular momentum loss rate is. k x 24/25 (eq. I29> . i n good 
agreemen t with the fo rce-free results found bv iGruzinovl J2005h 
and Spit kovskvl J2005h . 

(vi) For average integrated magnetization parameters as high as 
(To ~ 20, the radial position of the Y point (ry), which bounds 
the closed zone in the dipole models, is significantly less than the 
Light Cylinder distance. This result obtains despite the fact that the 
Alfven point rapidly approaches i?^ for cr > 1. The ratio ry/Ri^ 
is a very slowly increasing function of the surface magnetic field 
strength. An extrapolation of our results indicates that ctq must be 
very large, ctq s> 100, for the Y point to reach the Light Cylinder. 
That ry is generically much less than Rt in our calculations is due 
in part to the fact that although cto ~ 20 globally (Case C, ^5.3> the 



6.1 Rotation-Powered Pulsars 

Our results have several possible implications for rotation powered 
pulsars, where ctq is high. The fact that the Y-point is interior to 
the light cylinder (and the Alfven radius) suggests that observed 
braking indices less than 3 might be a consequence of equilibrium 
magnetospheric structure. We emphasize that the position of the 
Y point depends on a local equilibrium between the surface of the 
closed zone and the wind region. As a consequence, the value of the 
mass flux along the open field lines closest to the closed zone di- 
rectly affects the position of the Y point. The mass flux, in turn, de- 
pends on the boundary conditions at the stellar surface. In our sim- 
ulations, the mass loss rate is derived self-consistently, subject to a 
finite thermal pressure at r^vs ■ In the case of a pulsar, the injection 
of matt er at r^^ is thought to be due to non-MHD pair creation pro- 
cesses I Hibsch man & Aron s 2001) in a "gap" just above the stel- 
lar surface which, in an otherwise MHD flow, act to inject plasma 
with velocity already exceeding the sound speed. Then the mass 
loss rate is determined by the gap physics, not by the requirement 
of making a smooth transition from subsonic to supersonic flow. 
Thus, in principle ry depends upon the injection law determined by 
the pair creation physics at the surface. For this reason, we caution 
against over-interpreting our results |(vi)| and |(vTr)| in the pulsar con- 
text. Indeed, thermal pressure of the magnitude employed here is 
not likely to be relevant for classical pulsars.* On the other hand, 
magnetic dissipation and reconnection at and near the Y point might 
cause pressure and inertial forces to be significant in this localized 
region. In any case, our primary conclusions, that ry/R^ is gener- 
ally less than unity and that ry/Ri decreases as fl decreases, are 
intriguing possibilities now open to investigation with the advent 
of dynamical, large cr models of neutron star magnetopsheres. A 
more general study of the effect of the injection conditions on the 
structure of the magnetosphere and the accompanying wind is un- 
der way.' 

The appearance of outwardly propagating plasmoids at and 
beyond the Y-point as a consequence of (numerical) magnetic dissi- 
pation raises the intriguing possibility that noise in pulsar spindown 



* If the surface magnetic field is not highly stressed, the plasma in the 
closed zone is likely to be non-neutral and it will fill the magneto sphere via 
cross-field transport driven by shear flow turbulence i .Spitkovskv & AronsI 
2002), not via pressure, which alters the force balance at the Y point. 
' We are aware of the recent work of Komissarov (2005) who obtains 
ry « Rl for (T » 100. Because the injection conditions used by Komis- 
sarov (2005) are qualitatively different from the self-consistent calculation 
of the mass loss rate as a function of latitude obtained here, we do not be- 
lieve our two results are mutually contradictory. Instead, our finding that 
ry ai Rl only when (T 3> 100 serves to emphasize the point that the injec- 
tion conditions are critical in determining ry. 
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iCordes &Helf and"l980^ might arise from instability of the mag- 
netospheric currents due to real magnetic dissipation (e.g., via the 
coUisionless tearing instability). Our results show the possibility of 
20% or more torque fluctuations that could in principle give rise 
to a random walk in the rotation frequencies of pulsars, as is ob- 
served. Likewise, such torque fluctuations might give rise to noise 
and limit cycles in the observed phase of pulsars' subpulses, as is 
seen in many systems (Rankin 1986, Deshoandhe & Rankin 1999). 
Determining whether these observed phenomena could be due to 
magnetospheric dissipation requires treating the dissipation with a 
consistent physical model, which is an investigation beyond that 
reported in this paper. 

6.2 Proto-Neutron Stars & Proto-Magnetars 

High thermal pressure at the neutron star surface is thought to be 
a generic feature of neutron star birth and so the calculations pre- 
sented in this paper are directly relevant to young neutron stars, par- 
ticularly in the hypothesized rapidly-rotating and highly-magnetic 
initial state of magnetars. 

We identify five separate phases in any very young neu- 
tron star's life: (1) a pressure-dominated essentially non-magnetic 
phase in which the wi nd is driven by neutrino-heating (as in e.g., 
lOian & Wooslevll996h . (2) a phase in which magnetic field effects 
are present, but not dominant so that r^oj < 0. Ic « cy, where is 
the isothermal sound speed at the proto-neutron star surface, (3) a 
non-relativistic magnetically-dominated phase when Ra is greater 
than rjvs , but less than R^, (4) a relativistic phase in which ~ Rl, 
but ry < Rl, and lastly (5) an epoch when the force-free limit is 
applicable and ry - Ra - Rl- Phases (l)-(5) represent a time 
evolution starting immediately after the supernova explosion com- 
mences. The timescale for the evolution from phase (1) to phase (4) 
is set by the Kelvin-Helmholtz timescale for cooling of the proto- 
neutron star, Tkh ~ 10-100 seconds. The transition from phase (4) 
to phase (5) may occur on a longer timescale, or not at all, depend- 
ing upon the applicability of our results to classical pulsars as the 
MHD approximation breaks down. 

For parameters appropriate to a proto-magnetar, phase (3) lasts 
of order tkh- Our simulations show that in this phase, the wind is 
energetic and that the energetic flux is highly directed along the 
axis of rotation. The characteristic rotational energy loss rate in this 
phase is £ a 4 X lO'" Sj^Pf"*'^ ergs s"', where fin = 5(rjv5 )/10'* 
G is the "equivalent" monopole field (eq. l28> . and Pi = P/l ms. On 
the timescale Tkh, the total amount of energy extracted is compara- 
ble to the asymptotic supernova energy, ~ 10^' ergs. The magnitude 
of the rotational energy extracted in this phase and its coUimation 
along the rotational axis should have profound implications both 
for the spindown of millisecond magnetars and for the supernova 
remnants that accompany their birth. We save a detailed discussion 
for a future paper. 
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delay. These are corections due to the general relativisitc metric 
which are not present in the above cited papers {a(ri„) = 0.79). 

Notice that for the dipole, where the open field lines originate 
close to the pole, the base value of the magnetic field is about twice 
stronger than at the equator, and the rotational velocity is much 
smaller than the equatorial speed. 

This paper has been typeset from a TgX/ LSIgX file prepared by the 
author. 



APPENDIX A: PARAMETRIZATION OF THE MODELS 

All the simulations can be parametrized in terms of the ratios of 
characteristic velocities at th e injection radius (r , ,). Fol lowing the 
work by Mestel ( 1 968dh and lGoldreich & JulianI il970h . these are 
esentially the sound speed Ct, the rotational velocity v^, the non 
relativistic AL velocity -sJWJSnp, and the escape speed -\/2GM/r. 
For the 2D cases we consider the value of the rotational velocity 
and the magnetic field at the equator. Note that the characterisitc 
velocities defined above do not coincide with the cited works in 
newtonian gravity. For example the escape speed differs by a factor 
V2, and in the definition of one must include the effect of time 
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Table Al. Parametrization of the numerical models. 



ID Mod. 


AO 


A 


B 


C 


D 


E 


F 


G 


H 


1 


2GM 


1 1.44 


1 1.44 


6.74 


7.54 


8.98 


10.45 


1 1.44 


6.29 


6.29 


6.29 


b} 


061 


61 


71 


80 


95 


111 


60.6 


133 


533 


1330 


2GM 

r,-,. 


9.43 


9.43 


9.43 


9.43 


9.43 


9.43 


9.43 


9.43 


9.43 


9.43 


2D Monop. Mod. 


A 


B 


Bl 


C 


D 


E 










2GM 

2 

rifiCj 


1 1 AA 


1 1 AA 


1 1 AA 


1 1 AA 


1 1 AA 


1 1 AA 












0.61 


6.1 


6.1 


60.6 


606 


3030 










2GM 


9.43 


9.43 


4.20 


9.43 


9.43 


9.43 










2D Dip. Mod. 


A 


B 


Bl 


B2 


C 












2GM 


11.44 


11.44 


11.44 


11.44 


11.44 












If 


1.22 


12.2 


12.2 


12.2 


122 












IGM 


9.43 


9.43 


4.20 


38.5 


9.43 
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